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ABSTRACT 

We present the McMaster Unbiased Galaxy Simulations (MUGS), the first 9 galaxies 
of an unbiased selection ranging in total mass from 5x10^^ Mq to 2x10^^ M© simu- 
lated using n-body smoothed particle hydrodynamics (SPH) at high resolution. The 
simulations include a treatment of low temperature metal cooling, UV background 
radiation, star formation, and physically motivated stellar feedback. Mock images of 
the simulations show that the simulations lie within the observed range of relations 
such as that between color and magnitude and that between brightness and circular 
velocity (Tully-Fisher). The greatest discrepancy between the simulated galaxies and 
observed galaxies is the high concentration of material at the center of the galaxies as 
represented by the centrally peaked rotation curves and the high bulge-to-total ratios 
of the simulations determined both kinematically and photometrically. This central 
concentration represents the excess of low angular momentum material that long has 
plagued morphological studies of simulated galaxies and suggests that higher resolu- 
tions and a more accurate description of feedback will be required to simulate more 
realistic galaxies. Even with the excess central mass concentrations, the simulations 
suggest the important role merger history and halo spin play in the formation of disks. 
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1 INTRODUCTION 

Forming a galaxy like our own Milky Way remains a 
challenge for the currently accepted A Cold Dark Matter 
(ACDM) cosmogony. The Milky Way is comprised of three 
distinct, stellar components: a flattened, rotating disk] a 
compact, central and spheroidal bulge] and a diffuse, spher- 
ical halo of stars. Any consistent cosmogony needs to ex- 
plain the origin and evolution of each of these components. 
ACDM posits that the energy budget of the Universe is cur- 
rently dominated by vacuum energy (A), and the majority of 
the mass is invisible (dark) and only interacts with baryons 
via gravity. Thus, in the early Universe, thermal baryonic 
pressure did not support the dark matter, and because it 
is non-relativistic (cold) the dark matter first collapsed into 
small structures. Subsequently, the small structures merged 
hierarchically to form larger structures like the Milky Way. 

The ACDM paradigm provides a simple explanation 
for the formation of the stellar halo: stars formed early 
in small satellite galaxies, which got tidally stripped as 
their orbits brought them inside the tidal radius of the 
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main galaxy. While early observations indicated that stars 
in the Milky Way's halo might have condensed out of 
a monolithically collapsing gas cloud tegge n et aT]|l962l l. 
later observations found instead that formation through 
mergers like those propose d in the CDM paradigm are 
more likely (^Searle Zinni 19781 ). Today, digital surveys of 
the sky reveal structures in the Milky Way's stellar halo 
such a s streams and the remna nt cores of dwarf galaxies 
(iMaie wski 19931; | Belokurov et al . 2006). These are the exact 
signatures left by tidally disrupted satellites in simulations 
(Bullock Johnston 2005 : Abadi et al. 2006). 

Unlike halos, disks are not so neatly explained 
by ACDM. Although conservation of angular momen- 
tum naturally creates rotating disks, the hierarchi- 
cal buildup of structure impairs their formation, since 
disks form most efficiently in a quiet environment 
where gas cools and collapses smoothly. In ACDM- 
inspired simulations of substructure mergers, satellites or- 
biting disks ti dally he at sta rs turni ng thin disks into 
thick disks (iToth &; Ostrikerl Il99i lOuinn et al.1 Il993l : 



IVelazQuez k Whitelll999l : iKazantzidis et al .11 20081 ). Simula- 
tions of larger galactic mergers transform disks into centrally 
concentrated, spheroidal systems as the disks experience sig- 
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nificant angular momentum loss ( iBarnes &; Hernauis"tlll996l : 
ICox et all I2OO6I ). However, the observed distribution of 
galaxy morphologies can be reproduced with simple, ana- 
lytic models. In these models uniformly ro tating spheres col- 
lapse into centrifugally supported disks ([F^ 
198QI: iDalcanton et alJ ll997l : IMo et all 1 19981 : ' 



20061 : iGovernato et al.1 I2OQ7I: IScan nauieco et alJ ^ 200; 
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200 ih . The spheres start with angular momentum and mass 



profiles predicted using simulations of CDM structure for- 
mation. The contradiction between mergers and disk forma- 
tion may indicate that halo spin, not merger history, plays 
the dominant role in determining the morphology of galax- 
ies. 

The bulge is a spheroid of stars at the center of a galaxy. 
The spheroidal shape indicates that they formed as the re- 
sult of mergers. However, other evi dence revea l s tha t some 
bulges may have a secular origin. Kormendv (T^Fs) sug- 
gested that observations of rapidly rotating bulges indicates 
the existence of "pseudo-bulges", and recent simulations 
show that the central regions of isolated disks can buckle 
and cause stars to evolve th rough "peanut" shaped orbits 
into spheroidal distributions (Debattista et al]|2004l V 

Whether disks or spheroids form affects many galactic 
properties like their kinematics, color, light distribution, star 
formation history and metallicity in addition to morphology. 
Disk kinematics are dominated by ordered circular velocity, 
while in spheroidal components random velocity dominates 
over circular velocities. Because star formation usually hap- 
pens in galaxy disks where gas densities are sufficiently high, 
galaxies with more prominent disks have more recent star 
formation and display bluer colors, while spheroids tend to 
be more metal-rich and red. The prominence of morpholog- 
ical components also has an impact on the radial distribu- 
tion of galactic light profiles. Galaxies with more prominent 
spherical components exhibit more centrally concentrated 
light profiles. 

Models of galaxy formation require high resolution, hy- 
drodynamic numerical simulations. Analytic modeling can 
evolve a ACDM-motivated Gaussian density field into a 
spectrum of mass structures and populate those structures 
with stars such that they match the observe red lum inosity 
function of galaxies (I Cole et al ] [2OOOI : iBenson et al.1 l2003l : 
ISomerville et al.] I2OO8I I . However, because of the non- linear 
interactions of processes such as gas cooling, merging, tidal 
stripping, star formation, stellar feedback and active galactic 
nuclei, it is difficul t for these models to predic t gala xy mor- 
pho logies, though iBenson Devereuxl (|2009l ) and IPuttonl 
(|20 09) represent recent attempts. 

As an alternative, numerical simulations allow us to 
study how halos, disks, and bulges were created and evolve. 
Simulations that include gas are able to follow more physi- 
cal processes than simulations that only track the gravita- 
tional interaction of dark matter. Gas can be modeled using 
smoothed particle hydrodynamics (hereafter SPH), which 
partitions the gas in the Universe into particles and with a 
Lagrangian approach follows the motions of those particles. 
SPH is effective because it concentrates computational re- 
sources on the high density regions of a simulation, where 
galaxies form. 

Several studies of galaxies in a cosmological con- 
text have generated individual o bjects that are simi- 
lar to the observed lo c al galaxies (iGovernato et aP l2004l: 
[Robertson et al.1 l2004l : lOkamoto et al.1 l2005l : iBrook et all 
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I2OO9I : IPiontek & SteinmetJ 120091 : 
I2OO9I ). These require high resolution 



and large computational resources in order for several im- 
portant properties of the simulations to converge, such as 
the galactic structure, motions and star formation history. 
As a consequence of the high computational cost, the ini- 
tial conditions have been carefully chosen to maximize the 
chance that the simulation will produce the desired type of 
object (usually a disk galaxy). 

The previous cosmological simulations have shown that 
the collapse of gas is more complicated than smooth col- 
lapse into centrifugally supported disks. They show that a 
significa nt fraction of gas fiows into galaxies a long cold fil- 
aments ("Keres et al1l2005l : [Brooks et al.ll2009h . This raises 
the question: are centrifugal forces all that should support 
disks? Observations suggest that disks are supported by an 
equipartition of energy be tween ther mal, magnetic, and cos- 
mic ray pressure support (|Coxll2005l ). 



1.1 Stellar Feedback 

One way to introduce pressure support into simulations is by 
harnessing the energy massive stars release in stellar winds 
and supernovae explosions. Recent simulations continue to 
suffer from the overcooling that ha s long plagued morpho - 
logical studies of simulated galaxies. iNavarro Bend (|l99lh 
described how angular momentum transferred from gas in 
the disk to the dark matter causes excess central mass con- 
centrations, which leads to massive bulges that do not com- 
pare well to observations of disk galaxie s whose bulges ar e 
fainter and less massive than their disks ([Alien et al. I2OO6I ). 
"Navarro Bend ([l99l[ ) instead proposed that stellar feed- 
back could eliminate a significant amount of low angular 
momentum gas. They implemented a method to kinemati- 
cally excite gas particles around recently formed stars, but 
this showed littl e improvement in elimina ting low angular 
momentum gas ("Navarro Steinmetd [2"000i l. 

Stellar feedback plays a larger role in the development 
of satellite galaxies that merge with the parent galaxy. The 
effects of feedback in satellites can contribute to the final 
morphology of the main galaxy. If a satellite brings in more 
gas, it will contribute to make a larger disk; more stars will 
produce a larger spheroid. In the main galaxies, stellar feed- 
back may also determine the rate at which gas loses angular 
momentum and migrates through the disk into the dense, 
star forming center. 

Due to insufficient resolution, there is currently no sat- 
isfactory treatment for stellar feedback. Multiphase gas par- 
ticles attempt to capture the phenomenology of the IS M in- 
side individual particles ('Springel Hernauist[ [2003a[ l. but 
it is difficult to determine the appropriate pressure for each 
particle to exert on the others and s ometimes a stiff equa- 
tion of state needs to be enforced ([Springel &: Hernquist[ 
[2005[ ). If the gas dynamics are separated into hot and 



cold gas (Ritchie & Thomas[ [200l[ ). it is ambiguous when 
and how much gas should move from one phase to the 
other. Disabling the gas cooling reproduces s t ellar feedback 
([Gerritsen[[l997l : [Thacker Couchman[ [2OOOI : [Stinson et aP 
[2006[ 1. but resolution is often insufficient to turn off cool- 
ing in the proper amount of gas for the right length of 
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time, and SPH does not allow single particles to create 
their own outflow. Driven wi n ds reproduce galactic outflows 
(iNavarro Steinmetz 2000; Sprineel & Hernauist "2003b"; 
lODDenheime r & Dave 2006; Okamoto & Frenk 2009), but it 
is yet to be determined whether wind particles should be 
allowed to interact with surrounding gas and whether they 
provide suflicient pressure support to forming disks. Many 
potential avenues need to be followed to see how each feed- 
back recipe affects galaxy formation differently. To date, no 
stellar feedback model has been effective at reducing the cen- 
tral mass concentration in high mass galaxies, though strong 
adiabatic feedback combined with altering the star forma- 
tion threshold to a higher density in very high resolution 
simulations of low mass galaxies has had the m ost success 
(jMashchenko et aLllioOsI : iGovernato et al.ll2010h . 

Computational resources can now support surveys 
that include a range of galaxies simulated at moder- 
ate (1 million dark matter p a rticles inside Vyir) reso- 
lution (IScannapieco et aL I l2009l : lOkamoto FrenkI l2009l : 
IPiontek k, Steinmetj 12009 ) with different stellar feedback 
recipes. Each survey uses S PH, and each of t he previous 
surveys have used gadget. IScannapieco et al] (2009) sep- 
arated hot and cold gas and used super nova fee dback as a 



condu it from the cold to the hot ph ase. lOkamoto Frenk 
([2OO9I) used multiphase particles from ISpringel k, HernquistI 
(120 05') com bined w ith driven winds of different strengths. 
IScannapieco et al.l (|2009> ) found lowe r disk fractions than 
late ty pe spirals in their simulations. IPiontek k Steinmet j 
(|2009|) tested several stellar feedback recipes and did not 
find striking success with any of them. 

Observational samples of galaxies from galaxy redshift 
surveys, like the 2dFGRS (Colless et al. 2001) and SDSS 
(York et al. 2000), now contain millions of objects, allow- 
ing a much more complete view of not only the properties 
of typical galaxies, but also a quantification of how galax- 
ies are distributed within the multivariate parameter space 
of galaxy properties. In contrast, while many researchers 
are performing sophisticated galaxy formation simulations, 
each can only produce a handful of galaxies. Evaluating the 
success of these simulations requires a larger sample of sim- 
ulated galaxies. When only a small number of simulations 
exist, it is easy to find a good observational match for any 
one simulation; however, when a sample of simulated galax- 
ies exist that predict a mean and spread for any galactic 
property, these predictions can be directly confronted with 
the large observational samples. 

In order to address these problems, we have begun the 
McMaster Unbiased Galaxy Simulations (MUGS) project. 
The goal of MUGS is to generate a large sample of sophis- 
ticated galaxy formation simulations that sample potential 
sites of L* galaxy formation in an unbiased manner for direct 
comparison to the large observational samples now available. 
In this paper, we describe the methodology used for MUGS 
and present an overview of the first 9 simulations, particu- 
larly focusing on the relative formation of disks and bulges. 

MUGS provides an extende d look at galaxi e s sim - 
ulat ed using similar phys ics to iGove rnato et al.' (l2007l ) 
and IGovernato et al.] (|2009l l . Namely, supernovae are mod- 
eled with adiabati c "blastwave" feedback described in 
IStinson et al. M2006I ). In we describe how we created the 
initial conditions for MUGS, ^details the algorithm that 
evolves the simulations. SI examines the properties of the 



galaxies including their brightness, color, mass-to-light ra- 
tios, density profiles, bulge-to-total ratio, star formation his- 
tory, and metallicity. 



2 SIMULATIONS 

The simulations that comprise the MUGS sample each use 
the volume renormalization techniques from many previous 
simulations. The technique allows high resolution in a cos- 
mological context at reasonable computational cost. It fo- 
cuses resolution in one specific region of a cosmological vol- 
ume while simulating the rest of the volume at lower res- 
olution. The surrounding volume provides large scale den- 
sity waves and impart t idal torques on the region of interest 
(jOuin n Binnev 1992). 

We selected our galaxies from a cosmological cube 50 

Mpc on a side containing 256^ dark matter particles 
that was evolved to 2; = 0. The sim ulation is based on a 
4096^ realization of the CMBFAST (|Seliak k Zaldarriagal 
1996) power spectrum initially subsampled to 256^ to 
create a uniform resolution, dark matter-only volume. It 
uses a WMAP3 ACDM cosmology with = 73 km 
s~^ Mpc~\ Qrn=0. 24, Oa=0.76, Qbary=0.04, and (78=0.76 
(ISpergel et al.ll2007l ). 

The uniform volume was evolved to z—0 at which point 
the friends-of-friends algorithm was used to find the virial- 
ized halos with a linki ng length = inter-particle separa- 
tion (jDavis et al.l ll985V Every group with a mass between 
5 X 10^^ M© (705 particles in 256^) and 2 x 10^^ M© (2820 
particles in 256^) was examined to ensure that it did not 
evolve closer than 2.7 Mpc from any halo more massive than 
5 X 10^^ Mq. Massive structures contain hot gas that would 
significantly alter the evolution of a galaxy of interest. Simu- 
lating structures larger than 2 x 10^^ M© at the resolution of 
these galaxies is currently computationally unfeasible. Out 
of the 36,193 halos found with friends-of-friends, 761 were 
in the right mass range, and 276 of those were sufficiently 
isolated. From that sample, 9 halos were selected using a ran- 
dom number generator for more detailed simulation without 
regard for spin parameter or merger history. 

Particles within bvvir of each group's center at z=0 were 
traced to their positions in the initial conditions to specify 
the region of interest. In an effort to efficiently use computa- 
tional resources without loss of physical reality, the regions 
of interest have a non-spherical shape. Minimizing the num- 
ber of particles is critical since the n-body tree code must 
reconstruct the entire tree at each minor timestep, an op- 
eration that depends directly on the number of particles. 
The central region was filled with a regular grid of parti- 
cles to achieve an effective resolution of 2048^ at the cen- 
ter. Surrounding the non-spherical central region is a spher- 
ical region with a radius 1.2 times the maximum radius of 
the central region. This immediately surrounding region is 
populated with particles for an effective resolution of 512^. 
Outside this are three spherical regions equally spread in 
radius with effective resolutions of 256^, 128^, and 64^. The 
outskirts of the 50 Mpc cube are filled at an effective 
resolution of 32^. Each of these regions contains progres- 
sively more massive particles corresponding to the reduced 
resolution. Using this scheme, there are between 4 and 10 
times more particles in the high resolution region than in 
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ter only simulation. The values found when the galaxies are 
run at high resolution including baryonic physics vary mod- 
estly from these values. Figure [2] shows that in the initial 
sample of 9 galaxies does not include any from the high 
tail. As the sample grows in the future, it will more fully 
reproduce the multivariate distribution of halo properties. 
We note that only one halo with half mass redshift greater 
than 1.5 has a values greater than 0.05, but that some of 
the halos that acc rete ha lf of th eir mass later have A^ values 
that surpass 0.1. iRvdenI (| 19881 ) point out that late accre- 
tion contributes more angular momentum to halos because 
of their larger turn around radius. It is expected that galax- 
ies with the quietest merger histories will form galaxies with 
the most prominent disks. 




Figure 1. Projection of initial particle configuration for gl5784. 
The red particles represent gas and the light blue represents dark 
matter. Gas is only included in the high resolution region at the 
center. The dots become less dense further away from the center 
because of the decreased resolution and higher particle masses at 
larger radii. 



the surrounding low resolution regions combined. The reg- 
ular grids of particles in each region were perturbed using 
the Zel'dovich approximation with subsampled force resolu- 
tions matching the particle resolutions. This dark matter- 
only configuration was evolved to z=0. The particles that 
ended up inside Srvir in the dark matter-only resimulation 
had a fraction of their mass converted into a separate, neigh- 
boring gas particle with a mass corresponding to the cosmic 
baryon fraction, Qb/^m- We then reapplied the Zel'dovich 
perturbations to this new regular grid. Figure [1] shows the 
initial configuration for a sample galaxy. In the high resolu- 
tion region, dark matter particles have a mass of 1.1 x 10^ 
Mq and gas particles have an initial mass of 2.2 X 10^ M©. 
Stars form with a mass of 6.3 x 10^ M©. Each particle uses 
a gravitational softening length of 312.5 pc. 



2.1 Galaxy Sample 

One of the main goals of this study is to develop a bet- 
ter understanding of what sorts of galaxy morphologies are 
created using a wide range of initial conditions. The ran- 
domly selected MUGS sample spans a wide range of merger 
histories and halo angular momenta. The merger history is 
characterized as the redshift at which the galaxy obtained 
half of its final mass. The angular momentum of each halo is 
measured usins: A' = . (i Bullock et al.l l 200lh . Fig- 

^ ^y5/3GM3R ^ ^ 

ure[2] shows how the A^ and half mass redshifts of the halos in 
our 9 halo sample compare to the distribution of halos that 
passed our selection criteria. The A^ values shown in Figure 
[2] are the values from the low resolution, uniform dark mat- 



3 CODE 

The simulation s were evolved usi ng the parallel SPH 
code GASOLINE (jWadslev et al.1 1200^ 1. gasoline solves the 
equations of hydrodynamics, and includes radiative cool- 
ing. Gravity is cal culated for each particle using the 
(|l986l ) tree algorithm with tree elements that 
span at most — 0.7 of the size of the tree element's dis- 
tance from the particle, gasoline is multistepping so that 
each particle calculates its forces once every gravitational 
timestep Atgrav = ^\f^i where r\ — 0.175, ti is the grav- 
itational softening length (312.5 pc), and ai is the accel- 
eration. For gas particles, the timestep must also be less 
than Atgas = rycourant^, where ?7courant = 0.4, hi is the gas 
smoothing length and Ci is the sound speed. 

The cooling is calculated from the contributions of 
both primordial gas and metals as Atot p, T, Z) = 



:^{z, p,T). The primor- 



^0 

dial cooling follows the non-equilibrium evolution of inter- 
nal energy along with three primordial gas species (HI, 
Hel, and Hell). H2 cooling is not included. The scheme 
uses t he collisional ionisation rates reported in lAbel et al.l 
the radiative recombination rates from iBlackl (|l98lh 
and'Vern er &; Ferlandl (1 19961 ). and bremsstrahlung and line 
cooling from ICe nl (I1992I I. The metal cooling grid is con- 
structed using CLOUDY (version 07.02, last described 
bv Ferland et al. ( 1998 )), assuming ionisation equilibrium. 
A uniform ultraviolet ionizing backgrou nd, adopted from 
Haardt & Madau (in preperation; see iHaardt Madaul 
(1996)), is used in order to calculate the metal cooling rates 
self-consistently. The cooling lookup table is linearly inter- 
polated in three dimensions (i.e., p, z, T) and scaled linearly 
with metallicity. The energy integration independently uses 
a semi-implicit stiff integrator for each particle with the com- 
pressive heating and density (i.e. terms dependent on other 
particles) assumed to be constant over the timestep. 

The star formation and feedback recipes are the "blast- 
wave model" described in detail in IStinson et al.1 (I2OO6I ) 
with additional improvements as described in § 3.1. They 
are summarized as follows. Gas particles must be dense 
(rimin = O.lcm"^) and cool (Tmax = 15,000 K) to form stars. 
A subset of the particles that pass these criteria are ran- 
domly selected to form stars based on the commonly used 
star formation equation. 
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Figure 2. Three projections of the galaxy sample in mass, spin, and half mass redshift space. The light dots represent all the galaxies 
in our sample mass range from the uniform volume. The plus signs indicate the galaxies presented here. 



where is mass of stars created, c*" is a constant star for- 
mation efficiency factor, Mgas is the mass of gas creating the 
star, dt is how often star formation is calculated (1 Myr in 
all of the simulations described in this paper) and tdyn is the 
gas dynamical time. The constant parameter, ^ is tuned to 
0.05 so that the simu lated Isolate d Model Mi l ky W ay used in 
IStinson et aP (|2006l ) matches the Kennicutt (1998) Schmidt 
Law, and then c* is left fixed for all subsequent applications 
of the code. This star formation an d feedback treatment w as 
one of the keys to the success of iGovernato et aP (|2007l l in 
producing realistic spiral gala xies in acosmological simula- 
tion and the success of Brooks et aP (|2007l ) in matching the 
observed mass-metallicity relationship. 

At the resolution of these simulations, each star parti- 
cle represents a large group of stars (6.32 xlO^ M©). Thus, 
each particle has its stars partitioned into mass bins based on 
the initial mass function presented in iKroupa et al.l (|l993f ) . 
T hese masses ar e correlated to stellar lifetimes as described 
in iRaiteri et al.l (|l996l V Stars larger than 8 Mq explode 
as supernovae during the timestep that overlaps their stel- 
lar lifetime after their birth time. The explosion of these 
stars is t reated using the analytic m odel for blastwaves pre- 
sented in iMcKee Ostrikerl (|l977l ) as described in detail in 
IStinson et al. I (|2006l ). While the blast radius is calculated 
using the full energy output of the supernova, less than 
half of that energy is transferred to the surrounding ISM, 
EsN = 4 X 10^° ergs. The rest of the supernova energy is 
assumed to be radiated away. Iron and Oxygen are produced 
in SNII according to the analytic fits used in iRaiteri et al.l 
([199^: 

Mpe = 2.802 X 10""^ M^^-^^^ (2) 

Mo = 4.586 X 10"^M^^-^^^ (3) 

The iron, oxygen, and the supernova energy ejected 
from SNII are distributed to the same gas within the blast 
radi us. Each SNIa produces .63 M© Iron and 0.13 M© Oxy- 
gen ([Thielemann et al]|l986h and ejects it into the nearest 
gas particle for SNIa. 



3.1 Quantized Stellar Feedback 

One of the aspe cts of stellar feedback not given detailed 
consideration in IStinson et al L[2006) is the clustered na- 
ture of star formation. In Stinson et al. (2006), supernova 
feedback is chronologically distributed according to the stel- 
lar initial mass function and lifetimes. The combina t ion o f 
Padua group stellar lifetimes with the iKro upa et al.l (|l993h 
IMF results in a constant, small energy release each feed- 
back timestep (1 Myr, concurrent with star formation) for 
the 35 Myr until the lowest mass stars that explode as su- 
pernovae (8 Mq) explode. Since the blastwave radius and 
cooling shutoff time are calculated from the energy released 
during a timestep, the Stinson et al. (2006) method results 
in small blastwaves. iMcCrav k, KafatosI () 19871 ) describe how 
the accumulated energy of stellar winds and supernova feed- 
back create large superbubbles around star clusters. With a 
stochastic treatment of the energy release timing, we use the 
accumulated energy of all the stellar feedback that should 
result from a star particle to produce a larger, more realis- 
tic blastwave. The larger blastwaves should provide more 
pressure support to make more extended disks. Using a 
Kroupa et al. (1993) initial mass function, one supernova 
mass (> 8 M©) star is created for every 200 Mp of stars that 
form. Combining glob ular cluster (jllarrisll 19961) , open cluster 
( Piskunov et al.|[2008h , and embedded cluster (|Lada Ladal 
2003) catalogues, we estimate that a typical star forms in 
a cluster of 4000 M©. Thus, we require a minimum energy 
release of 20 supernovae (2 X 10^2 ergs) for the MUGS sim- 
ulations. 

We stochastically determine when a star particle re- 
leases feedback energy. The probability that energy is re- 
leased at each feedback timestep is 

^ _ {NsNii mod Nsnq) 

NsNQ 

where Nsnii is the number of supernovae calculated to ex- 
plode during that star formation timestep and Ns nq is the 
"supernova quantum", the number of supernova required 
per explosion. If the probability is greater than a random 
number selected between and 1, Nsnq supernovae's worth 
of energy is released. This causes SN energy to be released 
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sporadically over the 35 Myr until the largest star remaining 
is < 8 Mq. 



4 RESULTS 

In this section, we examine the galaxies that form from our 
set of cosmological initial conditions. First, we compare the 
properties of mock images of the simulated galaxies with ob- 
served galaxies. Second, we examine which parameters have 
the greatest impact on how much disk or spheroid forms. 
Third, we compare how much light is produced in our galax- 
ies with observations of the Tully-Fisher relationship and 
weak lensing mass measurements. Finally, we present the 
star formation histories and metallicity trends for the simu- 
lated galaxies. 

4.1 Tabulated Results 

In order to identify the main h alo and its satellites, we use d 
the Amiga Halo Finder (AHF) (iKnollmann k Knebdlioooh . 
AHF overcomes the difficulties Friends-of- Friends has in sep- 
arating neighboring groups, by instead identifying density 
peaks using an adaptive mesh algorithm. The adaptive mesh 
algorithm naturally leads to identification of substructure, 
which is critical for analyzing cosmological simulations. Ta- 
ble [1] summarizes key properties for each of the galaxies 
found with AHF in the sample. Galaxies are identified using 
the group number from the original friends-of-friends galaxy 
catalog. The columns are described below. 

The columns are defined as follows. Mass is the total 
mass in 10^^ M© located inside Vvir at 2; = of the final 
simulation including baryonic physics. is the spin param- 
eter of that matter defined as = . ^ zimm is 

the redshift of the last major merger, which is defined as 
the redshift when the AMIGA halo finder no longer distin- 
guished a satellite from the main halo where the satellite was 
greater than one-third the mass of the main halo at some 
point in its history. 2:1/2 is the redshift when the main halo 
was one-half its mass at ^ = 0. Mgas is the mass of gas in- 
side rvir at 2; = 0. /b is the baryon fraction of the final halo 
^ M^tar+Mgas \ ^^iQ mass of stars inside rvir at 2; = 0. 

Mdisk and Mtuige are the mass of stars kinematically classi- 
fied as part of the disk and bulge, respectively, as described 
in § 14.41 All the gas and stellar masses are reported in 10^° 
M0. The numbers of gas, star, and dark matter particles 
inside rvir at 2; = are Ngas, A^*, and Ndm, respectively. 

One of the problems inherent in running simulations 
where only a localized region is populated with high resolu- 
tion particles is that it is possible for low resolution particles 
or particles from outside the gas region to pollute the region 
of interest and cause unphysical results. No low resolution 
particles lie within rvir at 2; = of any of the MUGS sim- 
ulations. Some dark matter particles ("gasless DM") that 
originated outside the gas region ended up inside the virial 
radius. The presence of these particles without correspond- 
ing gas indicates that some of the gas in the simulation ex- 
perienced slightly less pressure than in reality. The initial 
conditions for gl536 were the first ones that were gener- 
ated, before the final criteria were firmly in place, and only 
contained gas particle pairs for dark matter particles that 



ended up inside rvir', as a consequence, it contains many 
more Bad DM particles than the other simulations. 

4.2 Simulated Images 

The most intuitive way to compare simulations with obser- 
vations is through mock images of the simulations. Such im- 
ages can be creat ed by assigning stellar population models 
like Starburst99 (Leitherer et al.l ll999f ) to each star parti- 
cle to determine the color and luminosity each star parti- 
cle should contribute to an image. Additionally, dust can 
modu late the image with extinction and scattering. Sun- 
rise (|jonssonl I2OQ6) is a Monte Carlo ray tracing program 
that produces simulated images by assuming dust exists in 
metal rich gas. Figure [3] shows images 50 kpc on a side that 
include scattering and absorption and are produced using 
SUNRISE. The image brigh tness and cont rast are scaled us- 
ing asinh as described in iLupton et al.l (2004) since disks 
have an exponential surface brightness profiles meaning the 
images span a wide range in surface brightness. 

Each galaxy is aligned so that the total angular momen- 
tum of the gas inside 1 kpc is pointed upwards in Figure 
[3l This presents the simulated galaxies edge-on to demon- 
strate the relative size of the disk and spherical component. 
In several of the images, a thin disk of young, blue stars is 
surrounded by a halo of old, red stars. In other images, lit- 
tle disk component is evident and the spherical component 
dominates. We note that several of the galaxies (g5664, for 
example) are not perfectly aligned indicating that the disks 
are warped or that the stars are orbitting around a slightly 
different axis than the gas. g7124 appears to be elongated 
vertically rather than horizontally because it is dominated 
by its spheroid which is strongly elongated perpendicular to 
its angular momentum of the inner gas, i.e. it is a spheroid 
rotating about the minor axis. 

In contrast to the images presented in Figure (3] the 
magnitudes and surface brightnesses used throughout the 
rest of the paper are derived from the face-on projection of 
the galaxies generated by sunrise, for which the extinction 
effects are minimized. 



4.3 Color-Magnitude Relationship 

One simple quantitative evaluation of the simulations is how 
the color and brightness of galaxies compare to observations. 
Figure |4] shows a g — r versus absolute r color-magnitude 
diagram (CMD) of - 3 X 10^ galaxies from the Sloan Dig- 
ital Sky Survey (S DSS) as the shaded two-dimensional his- 
togram (Bailin HarrisI boOS h . The colors and magnitudes 
have been inclination-corrected to their expected face-on 
values. The two well-known features of the observed CMD 
are the relatively narrow red sequence, which extends to 
very bright galaxies, and the more broad blue cloud, which 
is abruptly truncated at Mr ^ —23. Relatively few galaxies 
are observed to lie in the intermediate "green valley". The 
simulated galaxies are overplotted as the labelled points. 

The first conclusion that can be drawn from Figure |4] is 
that the simulated galaxies lie in observationally-populated 
regions of the CMD: they are representative of the colors and 
magnitudes of observed galaxies. Three simulated galaxies, 
g7124, g22795, and g25271, lie on the red sequence while the 
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Table 1. Simulation data 

Galaxy Mass 2;^^^ zi/2 h Mgas M* M^isk M^uige Ngas A^* Ndm Gasless 

(lO^i Mo) (all masses in lO^o M©) 10^ 10^ 10^ DM 



gl536 
g5664 
g7124 
gl5784 
g21647 
g22437 
g22795 
g24334 
g25271 



7.0 


0.025 


2.9 


1.8 


0.159 


5.1 


6.0 


1.8 


3.9 


2.4 


1.4 


5.3 


1640 


5.2 


0.025 


3.4 


1.3 


0.164 


3.8 


4.8 


1.3 


3.3 


1.8 


1.1 


4.0 





4.5 


0.039 





1.5 


0.163 


2.5 


4.8 


0.12 


3.2 


1.1 


1.1 


3.4 


147 


14 


0.0345 


3.42 


1.4 


0.150 


10 


11 


3.3 


5.5 


4.8 


2.6 


11 


2 


7.7 


0.048 


0.18 


0.53 


0.168 


5.8 


7.1 


1.1 


3.3 


2.7 


1.6 


5.8 


6 


8.8 


0.013 


1.56 


1.2 


0.170 


7.6 


7.3 


0.56 


5.6 


3.5 


1.6 


6.6 


1 


8.7 


0.011 


0.032 


1.3 


0.144 


6.0 


6.4 


0.26 


5.4 


2.7 


1.5 


6.6 


76 


7.7 


0.041 


0.085 


1.2 


0.162 


7.1 


11 


1.1 


5.4 


3.3 


2.4 


8.2 


60 


13 


0.014 


2.9 


1.3 


0.147 


8.9 


11 


2.2 


7.2 


4.0 


2.4 


10 


3 




-23 



-22 



-21 



-20 




-23 



-22 



-21 



-20 



Figure 4. Simulated galaxy g — r color as a function of absolute 
r magnitude in the Sloan Digital Sky Survey (SDSS) filter band- 
passes overlaid on an inclination corrected color ma gnitude dia- 
gram of SDSS galaxies from 2; < 0.2 from .Bailin Harris (^2008). 



Figure 5. As in Figure lU but with the SDSS sample restricted 
to isolated central galaxies with the same halo mass range as the 
simulated sample. 



remainder are members of the blue cloud. This is a slightly 
smaller proportion than the 48% of SDSS galaxies that lie 
on the red seque nce (within 0.08 of th e mode of the red 
sequence, i.e. the iBaihn k HarrisI (|2008h CMD parameter, 
CMD^ > -0.08). Given the small size of the simulated sam- 
ple, this difference should not be overinterpreted. However, 
a physical reason to expect different red sequence fractions 
is that the two samples are found in different environments: 
many observed galaxies are cluster members lying in mas- 
sive halos, while the simulated galaxies all lie at the centers 
of galaxy- mass halos. 

This is demonstra ted in Figure \E\ w here we have used 
the group catalogue of lYang et alj (I2OO7I I to restrict the ob- 
servational sample to central galaxies of halos with masses 
4 X 10^^ Mq ^ Mhaio ^ 2 X 10^^ Mq, and that do not have 
a neighbouring group with halo mass > 3 x 10^^ Mq within 
500 km s~^ and a projected radius of 2 Mpc. This is es- 
sentially identical to the selection criteria for our simulated 
galaxies. The relative fractions of blue versus red galaxies in 
Figure [5] are indistinguishable from those of our simulated 
sample: 30% for the observed sample and 3/9 = 33% for 
the simulated sample. It is also apparent from this Figure 
that not only do the simulated galaxies lie within the ob- 
served ranges of colour and magnitude for the global SDSS 
population, but they are also a particularly good match to 



observed galaxies that lie in halos of the correct mass and 
environment (although perhaps somewhat too luminous, a 
point we will return to in Section 4.7). 

While the MUGS galaxies are representative of some ob- 
served galaxies, they are unusually "green" : the red sequence 
galaxies lie on the blue edge of the sequence, while the blue 
cloud galaxies lie in the red half of the cloud. The simulated 
red sequence galaxies therefore have more star formation 
activity than usually observed, while the blue cloud galax- 
ies have either less recent or more ancient star formation 
than observed galaxies. One explanation for the sustained 
star formation in the red sequence galaxies is that it is a 
consequence of the environmental effects noted above. How- 
ever, the simulated galaxies are also too blue in Figure O 
which takes these environmental effects into account. Their 
blue color could also highlight some numerical failure of our 
simulations. MUGS does not include AGN feedback, which 
might be able to drive significantly more star forming gas 
out of the central regions of galaxies. The blue color may 
also result from overcooling, where lack of resolution causes 
excess gas to be driven to the galaxy center where it cools 
rapidly and forms stars due to the high gas density. 

The redness of the simulated blue cloud galaxies could 
also be a natural result expected for moderate mass galaxies 
evolving in an isolated environment, though Figure 5 again 
suggests that this is not the case. Alternatively, it could be 
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Figure 3. Edge-on mock images of the simulated galaxies using the g, r, and i filters found using the Monte Carlo radiative transfer 
program SUNRISE. Each image is a 2D projection of a box 50 kpc on a side. 



the result of excess ancient star formation due to our ne- 
glect of the increased feedback from metal-free Population 
III stars, or overcooling of gas that resulted in too many 
stars formed in the dense centers of galaxies at early times. 
Overcooling can in fact simultaneously make the blue galax- 
ies too red by building too large of a bulge (as seen in ^4.4|) , 
and make the red galaxies too blue by continuing to provide 
cold gas to the centers of bulge-dominated galaxies at late 
times. 

Regardless of these details, the colors and magnitudes 
of the simulated galaxies agree well with those of observed 



galaxies, giving us confidence that analyzing them will pro- 
vide a template of how galaxies form in the real universe. 

4.4 Bulge vs. Disk 

Another quantitative comparison between simulated galax- 
ies and observations is how the light and matter are dis- 
tributed. Figure [6] shows the face-on i band surface bright- 
ness profiles of the simulated galaxies each fit with the sum 
of a de Vaucouleurs r^^^ law (with effective radius, re) and 
exponential disk (with scale lengths, h) profiles. We calcu- 
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Figure 6. Surface brightness profiles of face-on mock i-band im- 
ages of the simulated galaxies, created using SUNRISE. The profiles 
are fit with the sum of an r^^^ law and an outer exponential disk. 
The bulge-to-total ratios are calculated using equation O 



late a ratio of the light from the bulge to the total light of 
the galaxy using 



B/T . 



(5) 



(|Binnev &: Merrifieldl Il998l ). The calculated B/T ratios 
shown in Figure [6] are all higher than 0.5, indicating that 
bulges dominate the emission from all of our galaxies. In 
contrast, observed B/T ratios for galaxies with classical 
bulge/disk profiles are often < 0.5, (see, e.g., figure 8 of 
Allen2006 base on the observed Millenium Galaxy Catalog) . 
Again, the high B/T ratios in the simulations indicate that 
too many stars form in the central region. 

Galaxies can also be separated into bulge and disk com- 
ponents based on their kinematics. For this analysis, disks 
are aligned such that the angular momentum of the gas 
within 30 kpc of the center of the halo points along the 
z-axis. Figure [71 shows the ratio of the 2;-component of the 
specific angular momentum vector to the total specific an- 
gular momentum for every star particle, as a function of 
its three dimensional radius from the center. Thus, Figure 
[71 represents a comparison between the angular momentum 
of the gas and the angular momentum of the stars. Each 
star within 40 kpc of the galactic center is classified as be- 
longing to the disk if greater than 0.75 of its total angular 
momentum is aligned with the gas, or to the bulge other- 
wise. The ratio of the masses of these components is given 
as the bulge to disk ratio. Note, though this is not the same 
bulge-to-total ratio derived photometrically as in Figure [6l 
the results are similar as the mass of stars in the bulge is 
always greater than the mass in the disk. Since the 40 kpc 
radius is much larger than what is generally considered part 
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Figure 7. The ratio of the z-component of the specific angular 
momentum vector to the total specific angular momentum for 
every star particle as a function of its three dimensional radius 
from the center. Stars are sorted into 100 x 100 bins 0.04 wide 
in log(r) and 0.02 high in jzr~^vj^. The horizontal dashed line 
at -^=0.75 indicates the separation between disk and spheroid. 

rvj_ ^ ^ 

The vertical dashed line at r=40 kpc indicates the radial extent 
inside which particles are classified as part of the spheroid. 



of the bulge, hereafter we refer to the kinematically defined 
component as the spheroid. 



4.4-1 What Creates the Disks? 

The bulge/disk decompositions show that MUGS simula- 
tions remain crude representations of the formation of disk 
galaxies; the bulge fractions are much larger than what are 
observed. Still, young stellar disks are apparent in at least 
gl536 and gl5784 in Figure [3l so we briefly examine which 
halo properties correlate with disk formation in the simula- 
tions. Figure [H shows the photometric bulge-to-total ratio for 
the galaxies as a function of mass, flnal spin parameter (A^), 
and the redshift of the last major merger. The use of unbi- 
ased simulations allows us to draw conclusions from these 
plots because no selection criteria was based on those vari- 
ables. One exception to our use of the photometric B/T ratio 
is g5664, which is the only galaxy with significantly different 
photometric and kinematic B/T ratios. The photometric de- 
composition for g5664 is unusual, and seems to display two 
distinct exponential components. The total bulge+disk fit 
is therefore poor and the decomposition based on that fit is 
suspect. Thus, for g5664, we have plotted the mean of the 
kinematic and photometric B/T ratios in Figure [8l 

A correlation is apparent between zimm and B/T. 
Galaxies with the most recent last major mergers have the 
largest bulges and typically the reddest colors in Figure [H 
One galaxy that had its last major merger long ago, but still 
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has a large bulge is g25271. We note that g25271 has one 
of the lowest spin parameter values. If we look in detail at 
the history of g25271, we find that it was involved in enough 
retrograde minor mergers that it never had a chance to grow 
a significant disk. 

We note that the galaxies with the most significant 
bulge separate into two groups, one with higher spin and 
one with lower spin. The higher spin galaxies all experienced 
recent last major mergers (also see D'Onghia & Navarro 
I2OO7I ). This indicates that it is only possible for galactic 
halos to gain angular momentum sufficient to increase its 
spin parameter above 0.04 if it has a late major merger. 
We further note that that these galaxies with a high spin 
parameter display the largest bulge fractions. Observed gi- 
ant low surface b rig;htness galaxies also exhibit large bulges 
(iLein et al.ll2010h . This might confirm the predictions of ana- 
lytic models of galaxy formation (e.g. Dalcanton et al. 1997) 
that halos with the highest spin parameters s hould host low 
surface brightness galaxies. iGovernato et al.] ((2009 ) showed 
how disks can form after late major mergers, but refrained 
from making any direct comparisons to LSBs. 

The lower spin group experienced a wide range of 
merger histories. Based on their spin parameters, it is ap- 
parent that they must have undergone a significant number 
of retrograde mergers. 

4.4-2 What Creates the Spheroids? 

Every galaxy in our sample displays a dominant spheroidal 
stellar component in both photometric and kinematic de- 
compositions. Since spheroids easily form as the result of 
merging, we examine the formation of the spheroid in the 
simulation that has the quietest merger history and largest 
disk (gl5784) here to discover whether there are any signifi- 
cant secondary effects. We leave study of spheroid formation 
in the other galaxies for future work. Figure [9] shows the or- 
bits of stars classified as spheroid stars in gl5784 using the 
kinematic decomposition. The stars are colored by age so 
that recently formed stars are yellow and stars that formed 
long ago are blue. The edge-on view shows that many of the 
recently formed "spheroid" stars are actually orbitting in 
the disk plane. The face-on view shows that many new stars 
are being formed in the central region of the galaxy with ra- 
dial velocities that lead to their classification as "spheroid" 
stars. While Figure [9] is colored to emphasize the young stars 
orbitting in a disk, the overall shape of the stars classified 
as part of the spheroid is spherical. 

Figure [10] shows the formation history of the stars in 
the spheroid at z=0. The spheroid stars are separated into 
three categories based on where they formed, 1) in satellites; 
2) within 1 kpc of the center of the main galaxy; 3) or in 
the disk of the main galaxy. We emphasize that stars in the 
disk category are no longer part of the disk, but are spheroid 
stars that formed in the disk. The decomposition was based 
on the distance stars formed from the center of the main 
halo. The stars were sorted into bins one million years in 
length based on their formation time. The center of mass 
of the stars that formed during that time was calculated 
by combining the stars formation locations in each time bin 
with the ce nter of the main halo found using the AMIGA 
halo finder ([Knollmann &: Knebei r2009) at the nearest out- 
put (generally full outputs were only written every 100 Myr). 




Figure 9. Velocities of stars kinematically classified as part of 
the spheroid in gl5784 at z=0 are shown as the thin lines. The 
stars are colored by their age with more recently formed stars 
yellow and less recently formed stars blue. The image is 80 kpc 
across. 
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Figure 10. Formation history of stars kinematically classified 
as part of the spheroid in gl5784 at 2; = 0. The solid line shows 
the total star formation history for all the stars that comprise 
the spheroid at z = 0. The blue dotted line shows the portion of 
that total that formed in the central 1 kpc of the main halo. The 
green dashed line shows the portion that formed in the disk and 
was heated to the spheroid. The red dot-dashed line shows the 
portion of the inner 40 kpc spheroid that formed in satellites. 

The large number of stars that formed within the central 
kpc of the main halo made finding the star formation center 
straightforward for the MUGS galaxies. We calculated the 
three dimensional distance between the center and the stars 
that formed during each time interval, and placed them into 
radial bins 1 kpc wide. Because the stellar density drops 
strongly at the edge of the disk, the existence of an empty 
radial bin serves as a delineation between the main galaxy 
and satellites. We therefore classify all stars within this ra- 
dius (but outside of 1 kpc) as "disk-formed" stars, and all 
stars outside of this radius as "satellite-formed" stars. The 
total bulge stellar mass is 1.1x10^^ Mo, the mass formed 
inside 1 kpc is 2.9x10^° M©, the mass formed in the disk is 
3.3x10^° Mo, and the mass formed in satellites is 4.1x10 
Mq. 

The final merger that disturbs stars from the disk into 
the spheroid happens about 13 Gyr into the simulation. 
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Figure 8. The photometric bulge-to-total ratios shown in Figure [6] as a function of galaxy mass, spin, and last major merger redshift. 
Galaxies at the bottom of these plots exhibit the strongest disks although they still have B/T > 0.5. 



After that time, nearly all the stars that constitute the 
spheroid form in the central 1 kpc. These are the stars that 
appear yellow in Figure (9] Stars forming in the unresolved 
central kpc is not limited to the final Gyr, they are a con- 
stant source of spheroid stars throughout the evolution of 
gl5784. For a significant fraction of the history of gl5784, 
the mass of stars formed in the central kpc of the main halo 
is similar to the mass of stars formed in satellites. While 
the formation of stars in the central 1 kpc is re miniscent of 
observations of rapidly rotating pseudo bulges (i Kormendvl 
'1993), we cannot draw any firm conclusions since the inner 
1 kpc is unresolved in these simulations. 

Stars that form in the disk also contribute a large frac- 
tion of stars to the final spheroid, though a merger 6 Gyr into 
the simulation marks the end of the dominant contribution 
of disk stars to the spheroid. We leave a detailed examina- 
tion of how stars migrate from the disk to the spheroid for 
future work, but speculate that most of the migration is due 
to the tidal disruption caused by merging satellites given 
the reduction in disk contribution to the spheroid o nce the 
merging activity is complete. We note that Debattista et al.l 
(|2004l l showed that secular evolution accounted for some 
migration of stars from the disk to the spheroid in high res- 
olution simulations of isolated disks. 

Satellites are expected to contribute a large mass of 
stars to the spheroid since stars are tidally stripped into 
spheroidal orbits. In the case of gl5784's quiet merger his- 
tory, the mass of the spheroid formed in satellites is similar 
to that formed in the inner 1 kpc and the disk. Finally, 
we note that the star formation in satellites is also largely 
confined to the unresolved central kpc similar to the main 
galaxy. 

The large number of stars that form in t he ce ntral 
kpc and disk echo the findings of IZolotov et all (|2009l ) and 
indicate that the predominant spheroid component is due 
in large part to stars that form in the main galaxy. How- 
ever, satellites contribute a significant mass of stars to the 
spheroid and disrupt a large fraction of stars out of the disk, 
so these simulations do not exclude traditional spheroid cre- 
ation through satellite mergers and tidal stripping. 



4.5 Density Profiles 

From the simulations, one can find the three dimensional 
distribution of all the matter in the galaxy as opposed to 
just that which is visible. Figure [TTI shows the resulting ra- 
dial density profile. Particles are sorted by radius and placed 
into 1000 bins, each containing the same number of particles. 
The dark matter only profiles for every galaxy are shown as 
the gray solid lines. The dark matter and total density pro- 
files are similar for each simulated galaxy. The dark matter 
profiles follow a power law or slightly shallower from the 
virial radius into O.lr^ir before mowing shallower, similar 
to dark matter only runs ('Navar ro et al.| [ 1993; [ Moore et al. 
1998; Reed et al. 2005). The total profile continues a steady 
rise along the power law into 1 kpc (about 0.005 Vvir) 
due to the stellar density before flattening out to a slope less 
than inside 1 kpc. Stars dominate the matter profile in- 
side of 2 kpc. 



4.6 Rotation Curves 

Rotation curves provide a comparison between the density 
profile found in these simulations and observations. Figure 
[T2l shows the circular velocity, Vc as a function of radius, 
where Vc = \/~^^- Mock observations are not used because 
there is too much scatter in the star particle velocities at a 
given projected distance from the center of the galaxy. 

The excessive bulge illustrated in ^4.41 is again appar- 
ent in the rotation curves. Rather than exhibiting flat rota- 
tion curves as are observed, the simulated rotation curves all 
have a peak at the center due to a large central concentra- 
tion of mass. iGovernato et al.l (|2007f ) showed that increased 
resolution without any changes to star formation can also 
limit the central concentration of matter as the lessened im- 
pact of two body interactions with halo dark matt er parti- 
cles minimizes angular momentum losses. Governato et al.l 
(|2010i ) show how higher resolution and consequent modifi- 
cations to the star formation threshold density can limit the 
mass concentration in simulations of slightly less massive 
galaxies. 
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Figure 11. Particles are sorted by radius and placed into 1000 
bins each containing the same number of particles. The radius 
of the bin is plotted as the mass weighted mean of the particle 
radii for each bin. The total density profile (various colours) for 
each halo is compared with the profile of the dark matter from 
the same simulations (thin grey). The arrow is placed at a radius 
of 2esoft^ the radius inside which IPower et al] (|2003l ) determines 
density profiles can no longer be trusted. 
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Figure 12. Rotation speed as a function of radius. The circular 
velocity (solid) defined as Vc = '\f^~ in radial bins for each of 
our galaxies. 



4.7 Tully-Fisher 

The Tully-Fisher (TF) relationship shown in Figure [13] com- 
pares the luminosity of a galaxy with its rotation velocity. 
Velocity indicates the dynamical mass for the inner parts 
of a galaxy, so the TF relationship indicates how bright a 
galaxy of a given mass should be. Most galaxies in the range 
of masses that we simulated have flat rotation curves, so the 
radius at which the rotation velocity is measured is irrele- 
vant. In the simulations, however, the rotation curve rises 
sharply initially and t he declines gradually w ith radius. We 
follow the example of iGovernato et aP (|2007ll and me asure 
the rotation velocity at 3.5/i. IGovernato et al.l (|2007f ) find 
in a resolution study that this is the radius at which ve- 
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Figure 13. I band absolute magnitude as a function of circu- 
lar velocity at 3.5 disk scale lengths (h). The large dots repre- 
sent the simulated galaxies while the small dots are from the 
iGiovanelli et al. | (1199711 sam ple of cluster galaxies and the trian- 
gles are from Icourteau et^al.. (,2007. ) . 



locity converges. It is unclear from Figure [12] that the ro- 
tation curves are flat at this radius, but they are flatter 
than they are at smaller radii. The I band magnitude for 
the galaxies generally matches the observations of the cor- 
respo nding; velocities, t houg jh the velocities remain slightly 
high. IGovernato et al.l (|2010l ) showed that higher resolution 
and higher star formation density threshold produce galaxies 
with more slowly rising rotation curves, so higher resolution 
simulations of these galaxies is an obvious next step. 



4.8 Mass-to-Light Ratio 

The relationship between the luminosity of a galaxy and the 
mass of its halo is a key observational question that provides 
constraints on semi-analytic models and strongly affects the 
interpretation of observational samples with respect to theo- 
retical predictions about halos. In Figure [Ml we have plotted 
the luminosity of our simulated galaxies versus their total 
(dark plus baryonic) mass. The luminosities are K-corrected 
to the °"^r band (the r band redshifted to z = 0.1) using 
KCORRECT v4.1.4 ( Blanton et al. 2003 ) using the SUNRISE- 
generated SDSS gri magnitudes in order to facilitate making 
comparison to the observations. 

By construction, our galaxies inhabit a relatively nar- 
row band in halo mass, from ^ 5 x 10^^ Mq to ^ 2 x 
10^^ Mq. Their luminosities cover a similar relative range 
in luminosity, from ^ 1.5 x 10^° Lq to 5 x 10^° Lq, for 
a typical total mass-to-light ratio of ^ 30 within the virial 
radius. 

We have overplotted a variety of observational estimates 
of the relationship between central galaxy luminosity and 
the total halo mass onto Figure [M] These observational tech- 
niques are: weak gravitational lensing (*Mandelbaum e t al.l 
2006, blue dot-dashed line for exponential-dominated late- 
type galaxies, and red dashed line for de Vaucouleurs- 
dominated early- type galaxies), a halo model simultane- 
ously fit to ga l axy c lustering and the luminosity function 
ICacciato^ et all (|2009l . solid green line), and the kinematics 
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of satellite galaxies iMore etal] (|2009l . dotted black line for 
the mean and shaded gray region for the measured spread 
in halo mass at a given luminosity). Note that while all of 
these studies are based on SDSS data, the techniques use 
completely independent properties of the galaxies to deter- 
mine the halo mass. We have adopted h = 0.73, as used in 
the simulations, to convert the /i- independent units used in 
these studies into physical units for direct comparison. 

When comparing the observations to the simulations, it 
is important to realize that in the simulations, the systems 
are selected based on the halo mass, while the luminosity is 
an output of the simulations; we have therefore plotted the 
halo mass as the independent variable and the luminosity 
as the dependent variable in Figure [M] However, in the ob- 
servational studies, the galaxies are selected based on their 
luminosity and the halo mass is what is measured; therefore, 
the luminosity is the independent variable and the halo mass 
is the dependent variable. This should be kept in mind when 
considering Figure 1141 and is why the weak lensing result 
around late- type galaxies appears mult i- valued. 

The halo mass-luminosity relationship for the simulated 
galaxies has the same shape and similar normalization as 
the observed relationship, pa rticu larly when comparing to 
the lMandelbaum et all (|2006l ) and lCacciato et all (|2009l ) re- 
sults. The simulated galaxies are, however, all more lumi- 
nous than the observed galaxies at a given halo mass, or 
equivalently the observed galaxies lie in more massive ha- 
los at a given luminosity. This difference is no larger than 
the difference seen between observational techniques, so it 
may not be significant. One aspect of the simulated and ob- 
served samples that may be important is that the simulated 
galaxies are constrained not to lie near large group or cluster 
halos; given the mass dependence of halo bias, those halos 
are themselves likely to be more massive. However, it seems 
unlikely that this is the explanation given th at the most 
isolated of the observational samples is that of iMore et al.l 
(^2009), with whom our results are the most discrepant. A 
more likely explanation is that this is another symptom of 
the overcooling problem, which turns too much of the gas 
in the system into stars and therefore results in an overlu- 
minous galaxy. 
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Figure 14. Galaxy luminosity in the ^-^r band plotted as a func- 
tion of the total halo mass including dark matter for the 9 simu- 
lated MUGS galaxies. Observational measurements of the mass- 
luminosity relationship are shown based on weak gravitational 
lensing (Mandelbaum et al. 200^) , halo model fits to th e galaxy 
clustering and luminosity function (ICacciato et al."2009V and the 
kinematics of satellite galaxies (jMore et al the measured 

dispersion is shown as the shaded region). 




4.9 Star Formation Histories 

One of the major ways to track the evolution of a galaxy 
is by examining when its stars formed. Figure [15] shows 
the star formation histories (hereafter, SFH) for each of 
our simulated galaxies. Peaks in the SFH correspond to 
merger events which drive starbursts. Mergers tidally dis- 
rupt disk gas and excite instabilities that cause gas o verden- 
sities (Mihos & Hernquist 1996; Springel 2000; Cox et al.1 
l2006l ). Following the merger peaks, the shape of the star 
formation history is an exponential increase followed by an 
exponential decline as the reservoir of gas is exhausted. The 
galaxies in Figure [15] that have the most and latest mergers 
correspond to those with the least well defined disks. 



4.10 Metallicity Evolution 

Trends in metallicity help track the star formation history 
of galaxies. Specifically, Figure [16] shows the distribution of 



Figure 15. The star formation rate plotted as a function of time 
for all the stars inside ryir (solid) and just stars that end up in 
the kinematically-defined disk (dashed). Stars are sorted by their 
formation time into 50 Myr bins and the total mass of stars in 
that bin is divided by 50 Myr. 



stars in [0/Fe] as a function of [Fe/H]. Evolution of metal- 
licity proceeds from low [Fe/H] as stars form and produce 
iron that is ejected during supernova explosions. Type II su- 
pernovae (SNII) occur on a shorter timescale than Type la, 
and they eject alpha elements such as carbon, oxygen, and 
silicon as well as iron. We use oxygen to track the abundance 
of alpha elements in our simulations. After stellar popula- 
tions age for 40 Myr, Type la supernovae (SNIa) start to 
explode, ejecting few alpha elements, but large quantities of 
iron. The range of [Fe/H] over which the [0/Fe] ratio re- 
mains high and constant indicates how many SNII explode 
before SNIa start contributing significant quantities of iron. 
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Figure 16. Metallicity distribution, [O/Fe] as a function of 
[Fe/H], of all the stars inside Tyir of the galaxy. Stars are sorted 
into 100 X 100 bins, 0.03 wide in [Fe/H] and 0.01 high in [O/Fe]. 
Evolution progresses from low [Fe/H] and high [O/Fe] to solar 
(0.0) [Fe/H] and [O/Fe]. The multiple evolutionary tracks through 
this diagram correspond to different stellar structures that each 
have their own star formation history. 



Thus, regions of active star formation like the Milky disk 
will generate long, high [O/Fe] plateaus befo re [O/Fe] de- 
creases due to the influ x of iron from SNIa (jBensbv et alj 
l2005l : iReddv et al.ll2006l ) . Places where star formation is not 
as efficient like dwarf galaxies exhibit lower [O/Fe] at low 
[Fe/H] distinct from abunda nces of stars measured in the 
Milky Way's disk and halo (|Venn et al.ll2004l ). Models of 
galaxy formation that couple N-body simulations with ana- 
lytic prescriptions for the stellar conte nt of satellite ga laxies 
are able to reproduce this dichotomy ( Font et al.ll2006l V 

Figure [16] shows that the majority of the stars in the 
simulations form with solar metallicity and abundances. At 
low [Fe/H], stars form with both high and low [O/Fe]. The 
peak [O/Fe] values are 0.4, which are lower than the ob- 
served abundan ces of t he di s k, which have been observed 
up to 0.6 Bensbv et al.1 (|2005l ): lReddv et alJ (|2006l l. This in- 
dicates a shortcoming of the simple power law fit used for 
oxygen enrichment. The fit only allows type II supernovae 
from stars up to 40 M© to produce oxygen, so the chemi- 
cal model does not capture chemical enrichment from more 
massive stars that produc e more Oxygen compared to Iron 
(|Wooslev k WeaverlfTooii ). 

Much like observational galactic archaeology, the dis- 
tinct metallicity evolution tracks apparent in Figure [16] pro- 
vide a useful alterative for discovering substructure inside 
galaxies as substructures exhibit stellar populations with dif- 
ferent metallicity signatures. We leave further discussion of 
such methods for future work. 



5 CONCLUSIONS 

We presented 9 galaxies from the McMaster Unbiased 
Galaxy Simulations simulated using N-body gravity and 
SPH. The galaxies are selected from a mass range around 
the mass of the Milky Way and from isolated environments, 
but their selection was otherwise unbiased for factors such 
as accretion history and angular momentum. 

The galaxies were examined using the radiative transfer 
program SUNRISE to enable comparisons between the simu- 
lated galaxies and real galaxies in the observed plane. The 
simulated galaxies have colors and magnitudes that compare 
well with a sample of inclination corrected isolated galaxies 
from SDSS, and in particular separate into the well-known 
red sequence and blue cloud. However, both simulated pop- 
ulations tend too much towards the "green valley" , indicat- 
ing that they contain more old stars than blue cloud galaxies 
and more young stars than galaxies along the red sequence. 

The surface brightness profiles of the simulated galaxies 
can all be fit with exponential disks combined with a central 
de Vaucouleurs r^^^ law similar to real galaxies. However, 
the proportion of bulge to total light (B/T) is higher than 
what is typically observed. The B/T ratios are also high 
when the stars are decomposed into the spheroid and disk 
based on their kinematics. There are no galaxies with a B/T 
fraction less than 0.5 whereas observed samples find many 
galaxies with B/T < 0.5. This result is similar to that found 
in many previous simulations. We note that many of the re- 
cently formed stars that are classified as part of the spheroid 
form with orbits in the disk plane in the central regions of 
the disk. We also note that most of the stars that comprise 
the spheroid formed in situ, but we leave the question of 
how the stars may migrate from the disk to a spherical dis- 
tribution for future work. 

As to the question of why galaxies form with more or 
less spheroid, there seems to be a modest trend with ac- 
cretion history. We find that the largest disks (gl536 and 
gl5784) form with a quiet merger history in which they 
had their last major merger prior to z—3, while the largest 
bulges all resulted from recent last major mergers. There 
also appears to be a dependence on halo spin as g25271 has 
a relatively quiet merger history, but low and results in 
a significant bulge and red color. Since all the galaxies used 
in MUGS fall in a limited mass range, these conclusions do 
not include the impact of mass. 

We compare the brightness of the final galaxies with 
their mass at two different radii in the final output. First, 
we compare our galaxies with the observed Tully-Fisher re- 
lationship that probes the amount of light a galaxy produces 
with the mass contained in its inner regions. Since the final 
mass concentration of all the galaxies is too high, we use 
a rotation velocity from 3.5 disk scale lengths away from 
the galaxy center. Previous studies have shown that this is 
the radius at which rotation velocities converge. The galax- 
ies are still slightly fainter than the observed sample based 
on these inner velocities. Second, we compare the bright- 
ness of our galaxies with observations of the whole halo 
mass derived using a number of different methods. In each 
case, the simulated galaxies are brighter than comparable 
observed galaxies at the same mass. This indicates that too 
many stars form in the simulation. The high central veloci- 
ties used in the Tully-Fisher relationship indicate that mass 
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gets too concentrated at the centers of the halos and while 
the galaxies are fainter than the observed TF relationship, 
the lack of resolution makes it difficult to determine whether 
the amount of stars formed is too many or too few. 

We are thus left with the challenge of creating more re- 
alistic simulations in order to obtain more accurate insights 
into how the important physical processes involved in galaxy 
formation result in the observed population of galaxies. For- 
tunately, there has been much recent work that guides the 
way forward. It has been sho wn repeatedly that higher res - 
olution makes better disks (jGovernato et al.l 1 20041 l2007l ) . 
More recently, it has been shown that high resolution com- 
bined with clustered star formation can r emove central den- 
sity concentrations from dw arf galaxies (|Mashchenko et aP 
I2OO8I : iGovernato et ahlboiOl l. 

While these simulations open many possibilities for ex- 
panding our understanding of how galaxies form, they also 
show that there is much work left to be done before we can 
claim to have simulated a sample of galaxies that compares 
well with real ones. 
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